Below is an intiial analysis of the DHW from the 2024 mass coral bleaching event on the Great Barrier Reef.
First, load the GBRMPA shp files for the Great Barrier Reef. For reproducibility the files are downloadable from Bozec et al 2022, and may differ from the official GBRMPA files available from the GBRMPA geoportal.
library(httr)
library(tidyverse)
library(janitor)
library(sf)
library(stars)
library(raster)
library(exactextractr)
library(tmap)
library(kableExtra)
library(reactable)
library(ggplot2)
library(terra)
library(tidyterra)
# write a function to get GBR shape file boundaries from eAtlas
get_gbr_shape <- function(crs="EPSG:20353", source="file"){
if (source=="url"){
#download from eAtlas
url <- "https://nextcloud.eatlas.org.au/s/xQ8neGxxCbgWGSd/download/TS_AIMS_NESP_Torres_Strait_Features_V1b_with_GBR_Features.zip"
destfile <- "TS_AIMS_NESP_Torres_Strait_Features_V1b_with_GBR_Features.zip"
GET(url, write_disk(destfile, overwrite = TRUE))
if (!dir.exists("unzipped_files")) {
dir.create("unzipped_files")
}
unzip("TS_AIMS_NESP_Torres_Strait_Features_V1b_with_GBR_Features.zip", exdir = "unzipped_files")
shapefile_path <- "unzipped_files/TS_AIMS_NESP_Torres_Strait_Features_V1b_with_GBR_Features.shp"
gbr_shape <- st_read(shapefile_path, quiet=TRUE) %>%
mutate(longitude = st_drop_geometry(.)$X_COORD,
latitude = st_drop_geometry(.)$Y_COORD) |>
filter(FEAT_NAME=="Reef") |>
st_set_crs(4283) |>
st_transform(20353) |>
st_make_valid() |>
clean_names() |>
dplyr::select(loc_name_s, qld_name, gbr_name, label_id, geometry, longitude, latitude) |>
mutate(gbr_name = as.factor(gbr_name)) |>
mutate(label_id = as.factor(label_id))
zones <- read.csv("/Users/rof011/GBR-coral-bleaching/data/gbr_zones.csv")
gbr_shape_output <- left_join(gbr_shape, zones, by="label_id")
return(gbr_shape)
} else if (source=="file"){
#read from disk
shapefile_path <- "/Users/rof011/GBR-coral-bleaching/data/Great_Barrier_Reef_Features.shp"
gbr_shape <- st_read(shapefile_path, quiet=TRUE) %>%
mutate(longitude = st_drop_geometry(.)$X_COORD,
latitude = st_drop_geometry(.)$Y_COORD) |>
filter(FEAT_NAME=="Reef") |>
st_set_crs(4283) |>
st_transform(20353) |>
st_make_valid() |>
clean_names() |>
dplyr::select(loc_name_s, qld_name, gbr_name, label_id, geometry, longitude, latitude,
area_ha) |>
mutate(gbr_name = as.factor(gbr_name)) |>
mutate(label_id = as.factor(label_id))
zones <- read.csv("/Users/rof011/GBR-coral-bleaching/data/gbr_zones.csv")
gbr_shape_output <- left_join(gbr_shape, zones, by="label_id")
return(gbr_shape)
}
}
### note that sf throws topology errors:
# Error in scan(text = lst[[length(lst)]], quiet = TRUE) :
# scan() expected 'a real', got 'TopologyException:'
# Error in (function (msg) :
# TopologyException: CoverageUnion cannot process incorrectly noded inputs.
gbr_shape <- get_gbr_shape(source="file") |>
mutate(id=sub("([a-zA-Z])$", "", label_id)) |> # drop multi-named sub-reef ID
group_by(gbr_name, id) |>
summarise() |>
st_make_valid()
# get distinct label_ID for each reef without sub-headings
# complex as some reefs (e.g. U/N) are catch-all "unknown reefs"
gbr_shape_label_ID <- get_gbr_shape() |>
as_tibble() |>
mutate(gbr_name=as.factor(gbr_name)) |>
mutate(id=sub("([a-zA-Z])$", "", label_id)) |>
dplyr::select(gbr_name, id) |>
distinct(gbr_name,id)
### this approach seems to give the same approx reef area as the raw data:
gbr_shape_raw <- get_gbr_shape()
# gbr_shape |> st_area() |> sum() # 28560112188 [m^2]
# gbr_shape_raw |> st_area() |> sum() # 28560118763 [m^2]
head(gbr_shape) |> kable("html") %>%
kable_styling("striped", font_size=11) %>%
scroll_box(width = "100%")
| gbr_name | id | geometry |
|---|---|---|
| (Big) Broadhurst Reef (No 1) | 18-100 | POLYGON ((1848551 7867941, … |
| (Big) Broadhurst Reef (No 2) | 18-100 | POLYGON ((1845249 7862813, … |
| (Big) Stevens Reef | 20-295 | POLYGON ((2094621 7653332, … |
| (Little) Broadhurst Reef | 18-106 | POLYGON ((1846710 7851185, … |
| (Little) Stevens Reef | 20-294 | POLYGON ((2085467 7656254, … |
| Abott Reef | 19-222 | POLYGON ((2119752 7720364, … |
There are 3724 reefs, with seperate polygons for each reef label
identifier (e.g. U/N Reef has 09-361a through to 09-361f). Individual
polygons per reef ID were merged with sf::union to form
multipolygons (see example for U/N 20-299 reef which has 9 separate
label_id components).
tmap_mode("plot")
a <- tm_shape(gbr_shape_raw |> filter(grepl("20-299", label_id))) +
tm_polygons("loc_name_s", legend.show=FALSE) +
tm_graticules() +
tm_layout(main.title="U/N Reef (20-299) ", main.title.position = "center", main.title.size = 1) +
tmap_options(max.categories = 57)
b <- tm_shape(gbr_shape |> filter(grepl("20-299", id))) +
tm_polygons("aquamarine4", alpha=0.2) +
tm_layout(main.title="U/N Reef (20-299) merged", main.title.position = "center", main.title.size = 1) +
tm_graticules()
tmap_arrange(a,b, ncol=2)

The NOAA Coral Reef Watch (CRW) daily global 5km satellite coral bleaching Degree Heating Week (DHW) is a metric that shows accumulated heat stress.
“There is a risk of coral bleaching when the DHW value reaches 4 °C-weeks. By the time the DHW value reaches 8 °C-weeks, reef-wide coral bleaching with mortality of heat-sensitive corals is likely. If the accumulated heat stress continues to build further and exceeds a DHW value of 12 °C-weeks, multi-species mortality becomes likely”
The data is available from 01/04/1984 to present day (with a ~2 day lag for processing). Data was extracted for cells within the GBR boundaries as determined by the boundary box (plus a 1km buffer) from the GBR shape file above. DHW were extracted for all 38 years (1985-2023) between 1st January to 1st July, on the assumption that the maximum heatstress would have passed by Austral winter (01/07) in each year. For each gridcell in each year the maximum DHW value was extracted to compare the DHW among years. DHW for 2024 were extracted seperately due to an incomplete timeseries at the time of analysis:
Note: each year of data is ~480mb in size, so total download size is ~15gb and may take some time
gbr_shape_bbox <- gbr_shape |> st_transform(4326) |> st_bbox()
gbr_shape_bbox <- gbr_shape |> st_transform(4326) |> st_bbox() |> st_as_sfc() |> st_buffer(1) |> st_bbox()
# custom function for DHW conversion to spatraster
# (there are better functions to download erdapp data but this works)
get_dhw_gbr <- function(timestart, timeend, timeout=60, crs="EPSG:20353"){
gbr_dhw_url <- paste0(
"https://coastwatch.pfeg.noaa.gov/erddap/griddap/NOAA_DHW.csv?", "CRW_DHW", "%5B(", timestart,
"T12:00:00Z):1:(", timeend, "T12:00:00Z)%5D%5B(",
-24.86075, "):1:(", -10.27304, ")%5D%5B(",
141.9369 , "):1:(", 153.2032 ,
")%5D"
)
temp_file <- tempfile(fileext = ".csv")
GET(url = gbr_dhw_url, write_disk(temp_file, overwrite = TRUE), timeout(timeout))
csvdata <- read.csv(temp_file, header = TRUE, skip=1)
write_csv(csvdata, paste0("dhw_", start_date, "_raw.csv"))
print(paste0("Downloading ", start_date, " - ", end_date))
csvdata <- csvdata |>
dplyr::select(2,3,4) |>
dplyr::group_by(degrees_east, degrees_north) |>
summarise(dhw = max(Celsius.weeks))
gbr_dhw <- tidyterra::as_spatraster(csvdata, crs="EPSG:4326") |> project("EPSG:20353")
return(gbr_dhw)
}
# Loop the function across years from 1981 to 2023
gbr_dhw_data <- list()
for(year in 1985:2023) {
start_date <- paste0(year, "-01-01") # Define start and end dates for the given year
end_date <- paste0(year, "-07-01")
dhw_data <- tryCatch({ # Fetch the DHW data for the GBR
get_dhw_gbr(start_date, end_date, timeout=500)
}, error = function(e) { #Use tryCatch to handle errors or timeouts
message(paste("Error fetching data for year", year, ":", e$message))
NULL # Return NULL if there was an error
})
gbr_dhw_data[[as.character(year)]] <- dhw_data #Store the list using the year as the name
}
# Ammend 2024 data
for(year in 2024) {
start_date <- paste0(year, "-01-01") # Define start and end dates for the given year
end_date <- paste0(year, "-07-01")
dhw_data <- tryCatch({ # Fetch the DHW data for the GBR
get_dhw_gbr(start_date, end_date, timeout=500)
}, error = function(e) { #Use tryCatch to handle errors or timeouts
message(paste("Error fetching data for year", year, ":", e$message))
NULL # Return NULL if there was an error
})
gbr_dhw_data[[as.character(year)]] <- dhw_data #Store the list using the year as the name
}
# function borrowed from FAMEFMR::saveSpatRasterList:
gbr_dhw_data_list<-rapply(gbr_dhw_data_list,f=terra::wrap, classes = "SpatRaster",how="replace")
qs::qsave(gbr_dhw_data_list,"gbr_dhw_data.RData")
As DHW is automatically calculated each day in the NOAA dataset, the
years in max DHW refer to the year of bleaching (e.g. 2016 = cumulative
heatstress from the end of 2015 and the start of 2016). For each reef
multipolygon in the gbr_shape dataset, the maximum DHW in
each year was calculated using an area-weighted approach due to
(exactextractr::exact_extract was used over
stars::st_extract | vect::extract |
raster::extract due to performance
issues]).
The exact_extract approach used below is conservative in
that for larger reef polygons (or multipolyons), it calculates the
mean value of cells that intersect the polygon weighted by
the fraction of the cell that is covered (as opposed to the
max value). Other metrics such as variance and
coefficient_of_variation may be of interest. Undefined (NA)
values are ignored in all of the named summary operations when they
occur in the value raster. When they occur in the weighting raster, they
cause the result of the summary operation to be NA.
library(reactable)
gbr_shape_dhw_list <- qs::qread("/Users/rof011/GBR-coral-bleaching/data/gbr_dhw_data.RData")
gbr_shape_dhw_list <- lapply(gbr_shape_dhw_list, function(x) terra::unwrap(x))
gbr_shape_dhw_list <- terra::rast(gbr_shape_dhw_list)
names(gbr_shape_dhw_list) <- paste0("year_", seq(1986,2024,1))
year_list <- paste0("year_", seq(1986,2024,1))
gbr_shape_dhw <- gbr_shape
for (i in year_list) {
year_col_name <- paste("dhw_max", sub("^.*_", "", i), sep="_")
gbr_shape_dhw[[year_col_name]] <- exact_extract(gbr_shape_dhw_list[[i]], progress = FALSE,
gbr_shape_dhw, fun = "mean", weights = "area")
}
head(gbr_shape_dhw) |>
kable("html") |>
kable_styling("striped", font_size=11) |>
scroll_box(width = "100%")
| gbr_name | id | geometry | dhw_max_1986 | dhw_max_1987 | dhw_max_1988 | dhw_max_1989 | dhw_max_1990 | dhw_max_1991 | dhw_max_1992 | dhw_max_1993 | dhw_max_1994 | dhw_max_1995 | dhw_max_1996 | dhw_max_1997 | dhw_max_1998 | dhw_max_1999 | dhw_max_2000 | dhw_max_2001 | dhw_max_2002 | dhw_max_2003 | dhw_max_2004 | dhw_max_2005 | dhw_max_2006 | dhw_max_2007 | dhw_max_2008 | dhw_max_2009 | dhw_max_2010 | dhw_max_2011 | dhw_max_2012 | dhw_max_2013 | dhw_max_2014 | dhw_max_2015 | dhw_max_2016 | dhw_max_2017 | dhw_max_2018 | dhw_max_2019 | dhw_max_2020 | dhw_max_2021 | dhw_max_2022 | dhw_max_2023 | dhw_max_2024 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| (Big) Broadhurst Reef (No 1) | 18-100 | POLYGON ((1848551 7867941, … | 0.4779795 | 6.827402 | 0.000000 | 0.1597545 | 0.7435788 | 0.0000000 | 0.000000 | 0 | 1.7756083 | 1.6970549 | 0.4418001 | 0 | 0.3921288 | 0.4099048 | 0 | 0.6378967 | 2.2812221 | 0 | 2.099609 | 0.6406564 | 0.0000000 | 0 | 0.0000000 | 0.0086903 | 0.6626197 | 0.2797973 | 0.8909236 | 0.9841175 | 0 | 1.9966893 | 2.639892 | 3.855308 | 0.0000000 | 0 | 4.845179 | 0 | 5.318794 | 1.0536289 | 3.319114 |
| (Big) Broadhurst Reef (No 2) | 18-100 | POLYGON ((1845249 7862813, … | 0.4809034 | 6.690235 | 0.000000 | 0.1626551 | 0.7322140 | 0.0000000 | 0.000000 | 0 | 1.8063999 | 1.7350723 | 0.3083107 | 0 | 0.2800516 | 0.4745746 | 0 | 0.5255354 | 1.9665437 | 0 | 2.207775 | 0.5768712 | 0.0000000 | 0 | 0.0000000 | 0.0397638 | 0.6499974 | 0.5694084 | 0.9334313 | 1.0662649 | 0 | 1.9251226 | 2.608260 | 3.985060 | 0.0000000 | 0 | 4.911906 | 0 | 5.216156 | 1.0033153 | 3.714746 |
| (Big) Stevens Reef | 20-295 | POLYGON ((2094621 7653332, … | 2.2126808 | 7.462231 | 0.000000 | 0.0000000 | 0.6870964 | 0.3540293 | 1.372475 | 0 | 0.0930394 | 1.1178683 | 0.0926058 | 0 | 1.7353671 | 0.1865064 | 0 | 0.0013572 | 7.0890689 | 0 | 4.832458 | 0.2291969 | 1.4887915 | 0 | 0.6001107 | 0.8347213 | 2.2275841 | 0.0000000 | 0.7364363 | 0.6408299 | 0 | 1.2891974 | 1.597199 | 4.002799 | 0.0894790 | 0 | 8.630629 | 0 | 5.692474 | 0.8342787 | 8.230654 |
| (Little) Broadhurst Reef | 18-106 | POLYGON ((1846710 7851185, … | 0.4927949 | 7.060939 | 0.000000 | 0.1851825 | 0.7358017 | 0.0000000 | 0.000000 | 0 | 1.7554580 | 1.9170403 | 0.6352469 | 0 | 1.0039911 | 0.5008817 | 0 | 0.5964878 | 2.9083798 | 0 | 2.061579 | 1.0029445 | 0.0133093 | 0 | 0.0000000 | 0.0000000 | 0.7555013 | 0.2112552 | 0.8949642 | 1.1583098 | 0 | 2.1319892 | 2.261860 | 4.152460 | 0.0000000 | 0 | 4.817797 | 0 | 5.554096 | 1.2604388 | 3.465622 |
| (Little) Stevens Reef | 20-294 | POLYGON ((2085467 7656254, … | 2.1869435 | 7.334218 | 0.000000 | 0.0000000 | 0.7461356 | 0.6261856 | 1.198692 | 0 | 0.0085639 | 0.9957811 | 0.1803312 | 0 | 1.9897678 | 0.2502062 | 0 | 0.0000000 | 6.7362270 | 0 | 4.782020 | 0.3121046 | 1.6231992 | 0 | 0.4489488 | 0.9118187 | 2.2421041 | 0.0000000 | 0.8632997 | 0.6447247 | 0 | 1.2262332 | 1.900885 | 3.752129 | 0.1179417 | 0 | 9.029800 | 0 | 6.001618 | 1.0063132 | 8.097749 |
| Abott Reef | 19-222 | POLYGON ((2119752 7720364, … | 0.7593529 | 5.319762 | 0.594146 | 0.1600000 | 0.3113776 | 0.0000000 | 0.000000 | 0 | 0.1686566 | 0.7345553 | 0.1500000 | 0 | 1.1293507 | 0.0000000 | 0 | 1.1302553 | 0.9077468 | 0 | 3.292808 | 0.1700677 | 0.4479577 | 0 | 0.3000000 | 2.9156933 | 0.7369853 | 0.4514952 | 0.4714453 | 0.0000000 | 0 | 0.4719019 | 2.083957 | 2.504015 | 0.0000000 | 0 | 4.460752 | 0 | 4.484212 | 0.6828330 | 9.708863 |
The output is an sfc of polygons for each reef with max
DHW for each year from 1984:2016 (see below example of the max DHW
experienced at Heron Island vs Lizard Island in 2016).
The DHW dataset can be filtered to quantify the trajectories of individual reefs through time (for example Heron Island in the southern GBR and Lizard Island in the northern GBR):
tmap_mode("plot")
a <- tm_shape(gbr_shape_dhw |> filter(grepl("Lizard", gbr_name)) |> mutate(dhw_max_2016=as.numeric(dhw_max_2016)) ) +
tm_polygons(col="dhw_max_2016", palette="-RdBu", breaks= seq(1:10), alpha=0.8, legend.show=FALSE) +
tm_layout(main.title="Lizard Island (2016)", main.title.position = "center", main.title.size = 1) +
tm_graticules(lwd=0.2)
b <- tm_shape(gbr_shape_dhw |> filter(gbr_name=="Heron Reef") |> mutate(dhw_max_2016=as.numeric(dhw_max_2016))) +
tm_polygons(col="dhw_max_2016", palette="-RdBu", breaks= seq(1:10), alpha=0.8) +
tm_layout(main.title="Heron Reef (2016)", main.title.position = "center", main.title.size = 1) +
tm_graticules(lwd=0.2)
tmap_arrange(b,a, ncol=2)

As the dataset is in sf format, the output for each year
can be plotted spatially to compare the maximum DHW for adjacent reefs,
for example Heron Reef and Lizard Reef (Lizard Reef is composed of
multiple polygons based on label_ID) where the within-reef trajectories
are fit with a simple GAM model:
dhw_subset <- gbr_shape_dhw |>
filter(gbr_name=="Lizard Island Reef (North East)" | gbr_name=="Heron Reef") |>
as.data.frame() |>
dplyr::select(-id, -geometry) |>
pivot_longer(-gbr_name, names_to="year", values_to="dhw") |>
mutate(year = as.numeric(str_extract_all(year, "\\d+"))) |>
arrange(gbr_name, year)
ggplot() + theme_bw() + facet_wrap(~gbr_name, ncol=2) +
geom_smooth(data=dhw_subset, aes(x=year, y=dhw, color=gbr_name),
method = "gam", formula = y ~ s(x, bs = "cr", k = 10), show.legend = FALSE) +
geom_line(data=dhw_subset, aes(x=year, y=dhw, color=gbr_name), show.legend=FALSE) +
geom_point(data=dhw_subset, aes(year, dhw, fill=gbr_name), shape=21, size=2, show.legend=FALSE) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
scale_fill_manual(values=c("aquamarine2", "coral2")) +
scale_color_manual(values=c("aquamarine4", "coral4"))

The dataset can be used to determine how the landscape has changed across multiple mass coral bleaching events since 1985 in terms of increasing frequency and intensity of heatwaves. The below example is the output of individual GAM model fits for 3724 reef polygons (label_ID) from 1985-2024.
# https://zenodo.org/records/5146061/preview/ymbozec/REEFMOD.6.4_GBR_HINDCAST_ANALYSES-v1.0.0.zip
# make a new zones for whole id
id_zones <- read.csv("/Users/rof011/GBR-coral-bleaching/data/gbr_zones.csv") |>
mutate(id=str_replace(label_id, "[a-zA-Z]$", "")) |>
dplyr::select(-label_id, -X) |>
distinct()
gbr_shape_dhw_centroid <- gbr_shape_dhw |>
st_centroid() |>
st_coordinates() |>
as.data.frame() |>
rename(lon=1, lat=2)
gbr_shape_dhw_zones_long <- gbr_shape_dhw %>%
as.data.frame() %>%
dplyr::select(-geometry) %>%
cbind(., gbr_shape_dhw_centroid) %>%
pivot_longer(cols = starts_with("dhw_max_"), names_to = "year", values_to = "dhw_max") %>%
mutate(year = str_replace(year, "dhw_max_", "")) %>%
left_join(id_zones, by="id")
ggplot() + theme_bw() +
geom_smooth(data=gbr_shape_dhw_zones_long |> na.omit() |> droplevels(),
aes(x=year, y=dhw_max, group=gbr_name),
color="darkgrey", se = FALSE, linewidth=0.05, method = "gam",
formula = y ~ s(x, bs = "cr", k = 10), show.legend = FALSE) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

Similarly, the data can be used to compare the impact of the current event (2024) with the previous mass bleaching events on the GBR (1998, 2002, 2012, 2016, 2017, 2020, 2022) in terms of maximum DHW experienced per event:
dhw_subset_bleaching_dhw <- gbr_shape_dhw |>
as.data.frame() |>
dplyr::select(-id, -geometry) |>
pivot_longer(-gbr_name, names_to="year", values_to="dhw_max") |>
mutate(year = as.numeric(str_extract_all(year, "\\d+"))) |>
arrange(gbr_name, year) |>
filter(year %in% c(1998,2002,2016,2017,2020,2022,2024)) |>
group_by(year) |>
mutate(mean_dhw_max = mean(dhw_max, na.rm=TRUE))
a <- ggplot() + theme_bw() +
ggridges::stat_density_ridges(data=dhw_subset_bleaching_dhw, alpha=0.8, scale=1.2, size=0.3,
aes(y=as.factor(year), x=dhw_max, fill=mean_dhw_max), rel_min_height=0.00001,
quantiles = c(0.5), quantile_lines = FALSE, quantile_fun = mean,
vline_color = c("white"), vline_width = 2, show.legend=FALSE) +
scale_fill_distiller(palette = "Reds", direction = 1, name="mean\nmax_DHW\nper reef") +
xlab("Maximum DHW per reef") + ylab("Mass bleaching event") +
theme(panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank())
a

# b <- ggplot() + theme_bw() +
# ggridges::stat_density_ridges(data=dhw_subset_bleaching_dhw, alpha=0.8, scale=0.8, size=0.3,
# rel_min_height=0.001,
# aes(y=as.factor(year), x=dhw_max, fill=mean_dhw_max),show.legend=TRUE) +
# scale_fill_distiller(palette = "Reds", direction = 1, name="mean\nmax_DHW\nper reef") +
# xlab("Maximum DHW per reef") + ylab("") + scale_x_continuous(limits=c(0,8), oob = scales::oob_squish) +
# theme(panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank())
#
# a+b+ patchwork::plot_layout(widths = c(0.7,0.3))
Previous studies have shown considerable variability in the spatial footprint of mass bleaching events on the GBR, with the southern region remaining relatively unimpacted by thermal stress in the past 4 mass bleaching events.
The sf output can be used to generate interactive maps
to explore the extent of DHW across reefs in the current 2024 event
(click individual reefs to see the reef_name / label_ID and DHW from
previous events in 1998, 2002, 2012, 2016, 2017, & 2020):
gbr_shape_dhw_map <- gbr_shape_dhw %>%
mutate(across(where(is.numeric), round, digits = 1))
tmap_mode("view")
tm_basemap("Esri.WorldImagery") +
tm_shape(gbr_shape_dhw_map) +
tm_polygons(col="dhw_max_2024", title="DHW Max 2024", palette="-RdBu",
breaks= seq(0,20,1), alpha=0.8, legend.show=FALSE,
popup.vars = c("dhw_max_1998", "dhw_max_2002", "dhw_max_2016",
"dhw_max_2017", "dhw_max_2020", "dhw_max_2022", "dhw_max_2024")) +
tm_view(set.view = c(145.3, -14.5, 10))